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A  FAST  METHOD  FOR  SOLVING 


A  CLASS  OF  TRI -DIAGONAL  LINEAR  SYSTEMS 

by 

Michael  A.  Malcolm  and  John  Palmer 


ABSTRACT 


The  solution  of  linear  systems  having  real,  symmetric,  diagonally 
dominant,  tridiagonal  coefficient  matrices  with  constant  diagonals  is 
considered.  It  is  proved  that  the  diagonals  of  the  LU  decomposition  of 
the  coefficient  matrix  rapidly  converge  to  full  floating-point  precision. 
It  is  also  proved  that  the  computed  LU  decomposition  converges  when 
floating-point  arithmetic  is  used  and  that  the  limits  of  the  LU  diagonals 
using  floating  point  are  roughly  within  machine  precision  of  the  limits 
using  real  arithmetic.  This  fact  is  exploited  to  reduce  the  number  of 
floating-point  operations  required  to  solve  a  linear  system  from  8n-7 
to  5n+2k-3  ,  where  k  is  much  less  them  n  ,  the  order  of  the  matrix. 

If  the  elements  of  the  sub-  and  superdiagonals  are  1  ,  then  only  4n+2k-3 
operations  are  needed.  The  entire  LU  decomposition  takes  k  words  of 
storage,  and  considerable  savings  in  array  subscripting  are  achieved. 

Upper  and  lower  bounds  on  k  are  obtained  in  terms  of  the  ratio  of  the 
coefficient  matrix  diagonal  constants  and  parameters  of  the  floating-point 
number  system. 

Various  generalizations  of  these  results  are  discussed. 


1.  Introduction 


We  will  consider  the  solution  of  linear  algebraic  systems  having 
real  symmetric,  diagonally  dominant,  tridiagonal  coefficient  matrices 
with  constant  diagonals.  This  problem  occurs  frequently  in  solving  certain 
kinds  of  partial  differential  equations,  boundary  value  problems  of  ordin¬ 
ary  differential  equations,  and  cubic  spline  interpolation  problems. 

Consider  the  coefficient  matrix 

m 

b  i 

a  b 
b  a  b 

•  .  •  > 

•  •  • 

b  a  b 
b  a 

m 

of  order  n  .  The  usual  LU  decomposition  of  A  requires  n-1  divisions, 
n-1  multiplications,  and  n-1  additions.  The  solution  of  the  equations 
LUx  =  d  requires  an  additional  n  divisions,  2n-2  multiplications,  and 
2n-2  additions.  With  the  following  observation,  the  entire  LU  decompos¬ 
ition  of  A  can  be  stored  in  k  floating-point  words,  and  the  solution 
of  the  linear  system  Ax  =  d  can  be  obtained  in  k  divisions,  in-l  mul¬ 
tiplications,  and  2n-2+k  additions,  where  k  is  usually  much  less  than 
n  .  Typically,  k  is  on  the  order  of  10.  Moreover,  k  can  easily  be 
estimated  from  the  values  of  a  and  b  and  parameters  of  the  floating¬ 
point  number  system  used  in  the  solution.  If  b  =  1  ,  then  n  multiplies 
can  be  avoided.  In  addition  to  a  smaller  operation  count,  substantial 
savings  in  array  indexing  are  achieved. 


^ .  The  Algorithm 

Consider  the  matrix 


a  1 

1  a  1 

1  a  1 

B  =  I  «  •  • 

•  *  4 
.  •  * 

1  cr  1 

1  a 


where  a  =  a/b  .  Note  that  A  =  bB  .  The  analysis,  as  well  as  the  comp¬ 
utation  is  simplified  by  considering  the  coefficient  matrix  to  be  B  and 
the  linear  system  bBx  =  d  .  B  can  be  factored  into  the  product  LU  , 
where 


L  = 


m  m 

- 

l 

u.  1 

1 

fL  1 

"2  \ 

• 

t2  1 

,  u  = 

• 

• 

% 

,  •  , 

* 

1 

• 

• 

4 

• 

l  ,  1 

u 

n-1 

n 

to  - 

m 

using  the  recurrence  relations; 

u^  =  o  >  =  Vui_i  >  u.  =  <y  - 


or 


u.  =  a  -  l/u.  ,  ,  i=2 , . . . ,n 

i  '  l-l  ’  ’ 


(1) 


2 


Under  suitable  conditions,  to  be  discussed,  the  converge  and 

=  i^+i  =  •  •  •  =  in  =  £  to  machine  accuracy.  In  the  computer,  one  simply 
computes  and  stores  the  values  of  t  ,  i»l,...,k  .  The  solution  vector 

x  can  then  be  computed  as  follows: 


*1  =  di  * 


y  =  d,  -  £i  1yi  1  ,  i=2,...,k  ,  y.^  =  d.^  -  iy^,  i=k+l,...,n  , 


zn  =  £yn  * 


zi  “  *(yi  “  zi+l^  »  i^-1*-***11  »  zi  =  ^i^yi  “  zi+l^5  1  > 


-1 


x.  -  b  z.  ,  i=l, . . . ,n  . 
i  X  '  * 


3 .  Convergence  of  the  LU  decomposition 

We  will  show  that  when  A  is  diagonally  dominant,  the  sequences 
[u.^]  and  [£ ^]  converge.  We  will  also  find  an  estimate  of  the  rate  of 
convergence  which  can  be  used  to  determine  a  value  for  k  . 

It  is  sufficient  to  show  that  the  sequence  [u.,]  converges,  and 
for  this  we  assume  diagonal  dominance,  o^  equivalently,  \a\  >  2  .  The 
following  theorem  is  a  special  case  of  a  theorem  of  Parter  (1962)  for 
band  matrices. 


Theorem  1:  If  | cr |  >2  ,  then  the  sequence  [u^]  converges  to  u  where 


a  +  sgn(g) 
2 


Proof;  Convergence  follows  from  the  fact  that  the  sequence  [aru^]  is  bounded 
and  monotone; 


(2) 
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Lemma  1  (boundness):  If  |<*|  >  2  ,  then 


a\ii  >  2  ,  i=l,.. 


Proof:  From  (l),  u1  =  or 
for  some  value  of  i  >  1  , 


Thus  oru^  =  a  >4 
By  (i), 


(4) 


Now  assume  that  (4)  holds 


cm 


2  2 


i+1 


=  a  -a/*  >  a  -  a  /2  >  2 


Lemma  1  follows  by  induction. 


Lemma  2.  (monotonicity):  If  a  >  2  ,  then 


cmi+1  <  cn^  ,  i+1, ...  . 


Proof:  From  (l), 


“2  *  “l  '  -  -  > 

and  «(uJ+1  -  Uj)  =  ■■■  «(Ul  -  *U1)  , 


i  i-1 


It  follows  from  Lemma  1  that  the  u^  must  all  have  the  Same  sign.  Thus, 


by  induction, 


“(ui+i  -  ui>  <  0  ■ 


Now,  in  the  limit, 


u  =  O'  - 


or. 


u  -  cm  +  1  =  0  . 


Equation  (.i)  is  the  quadratic  formula  with  the  sign  of  the  radical  chosen 
to  avoid  a  contradiction  with  Lemma  1.  This  completes  the  proof  of 
Theorem  1. 
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The  following  two  theorems  provide  a  way  to  estimate  the  value  of  k 


Theorem  2;  If 


k  < 


a  >  2  ,  then 


t  -  1  - 

logp(c£ 


(5) 


where  p  is  the  floating-point  radix,  t  is  the  number  of  digits, 
and  f«l  denotes  the  smallest  integer  not  less  than  t  . 

Proof :  We  will  first  prove  the  following  lemma. 


Lemma  3 :  If  | a |  >  2  ,  then 

2 


a(ui+l  ‘  V  >  ~  ~u~  ~  i)1"1*  i=1»-* 


Proof;  From  (l).  Lemmas  1  and  2, 

"("w  -  ui>  ?  u.u.  .  “(ui  •  Ui-1>  <  0  ’ 
1  1-1 


and 


uiVi  mi  -  1 


>  0  ,  i=2 , . . .  . 


2  a  i  o 

Now,  aru  =  a  — - -  ,  i=2,...  . 

1  ui-l 


By  Lemma  2,  and  the  fact  that  u.^u  >  0  , 


a  <  a  i=l 
11  '  * 


u.  u 

i 


Thus,  cm.  >  a2  — —  , 
1  u 


and 


“S - T 


^  -  -r  - 1 


,  i=*2 , . . . 


Thus,  a^u..,  -  u  )  > 


i+1  “i'  '  2  a  . 

or - 1 

u 


o?(ui  -  u^_^)  >  1=2 »' 


(6) 

(7) 

(8) 

(9) 
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Repeated  application  of  this  inequality  yields 


(Ui+L  -  \)  >  (<f  -  -  l)1-1a(u2  -  ux),  i=l, . 


Since  a(n^  -  u^)  =  -1  , 
the  Lemma  is  proved. 

Dividing  (6)  by  cm  >  0  and  talcing  absolute  values, 


(10) 


Ui+1  '  Ui 


u 


<  -i-  to2  -  JL  -  I)1*1 

cm  u 


(11) 


1-t 


Requiring  the  right-side  of  (ll)  to  be  less  than  p  "  gives  a  sufficient 
condition  on  i  for  the  convergence  of  [u^]  .  Taking  logarithms  yields 


the  sufficient  condition 

t  -  1  -  log  cm 


i  >  1  + 


l°6p(«2  -  -  1) 


(12) 


Thus  k  need  be  no  larger  than  the  smallest  possible  value  of  ji^  given 
by  (12). 

Theorem  3:  If  |or|  >  2  ,  then 

log  aru 

J _ 


t  -  1 

>  1  + - 


W- 2) 


Proof ;  We  will  first  prove  the  following  lemma. 


Lemma  4;  If  |  a  |  >  2  ,  then 

a(ui+l  "  Ui}  -  "(Qr2  '  2)1_i>  i=1”-*  *  (13) 

Proof :  By  Lemma  2  and  (l), 

.  2  <  , 
cm^  <  oru1  =  a  ,  i=l ,  —  . 
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» 


Since,  by  Lemma  1,  cm^  >  0 


a 


u.  — 
1 


>  1  ,  i=l, . . , 


Substituting  into  (8)  and  (9)  gives 


1  >  -  ,  i=2 , . 


uiui-l  a2  -  2 


This  inequality  and  (7)  and  (10)  yield  Lemma  1+. 

Dividing  (12)  by  am  >  0  and  taking  absolute  values  gives 


Ui+1  "  Ui 


u 


1  /  2  „  vl-i  .  . 

> -  (a  -  2)  ,  i=l,...  . 

-  cm  v  * 


(15) 


Setting  the  right-side  of  (13)  greater  than  gives  a  non-convergence 

condition  for  i  ,  and  thus,  a  lower  bound  on  k  .  Taking  logarithms 
yields  Theorem  3.  | 

If  we  denote  by  K  ,  the  upper  bound  given  in  Theorem  2,  and  by 
K  ,  the  lower  bound  given  in  Theorem  3,  we  have 
K  <  k  <  K  . 

In  practice,  these  bounds  are  very  close.  Usually  K  =  k  =  K  .  The 
following  table  gives  values  for  K  ,  K  and  k  for  various  values  of 
a  for  both  single  and  double  precision  on  the  IBM  360. 
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% 


or 

:  S 

X 

hort  Precis 

(p=l6,  t=6) 

k 

ion 

K 

f  1/311 

(f}= 

K 

g  Precisio 

16,  t=l4) 

k 

n 

K 

2.05 

18 

27 

30 

46 

77 

80 

2.1 

16 

20 

22 

4l 

55 

57 

2.2 

14 

15 

16 

35 

4o 

4l 

2.3 

12 

13 

13 

31 

33 

34 

2.4 

11 

11 

11 

28 

29 

29 

2.5 

10 

10 

10 

25 

26 

26 

3.0 

8 

8 

8 

19 

19 

19 

4.0 

6 

6 

6 

14 

l4 

14 

5.0 

5 

5 

5  ? 

12 

12 

12 

6.0 

4 

4 

4 

11 

11 

11 

7.° 

4 

4 

4 

10 

10 

10 

Upper  end  Lover  Bounds  (K  end  K)  end 
Observed  Values  for  k  for  the  IBM  360 
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The  preceeding  theorems  characterize  the  convergence  of  the 
sequence  [u^]  in  the  absence  of  rounding  errors.  If  the  computer  arith¬ 
metic  satisfies  certain  reasonable  rules,  then  the  computed  sequence  [u^] 
also  converges  monotonically  to  a  limit  u  which  is  very  close  to  u  . 

We  will  prove  this  result  for  a  >  2  .  A  similar  argument  holds  for 
a  <  -2  . 

Let  denote  the  operation  of  floating-point  divide,  and  © 
denote  the  operation  of  floating-point  subtraction.  For  any  floating¬ 
point  numbers  a,  b,  and  c,  we  will  assume  the  following: 

(i)  a  >  0  3  10a  >  0 

(ii)  a>b>l3l>  10b  >  10a 

(iii)  a  >  b  3  c@b  >  c@a 

(iv)  a  >  2  3  a01  >  1 

(v)  a@0  =  a 


Theorem  4 :  If  a  >  2  ,  and  the  computer  arithmetic  satisfies  the  above 
rules,  then  the  computed  sequence  [u^]  converges  monotonically  to  u  and 
u  =  u  +  0  (P1-t)  . 

Proof;  u^  =  a  >  2  and  =  "-0(1©  a)  .  Since  a  >  2  ,  (i)  yields 
10«>O.  From  (iii)  and  (v)  we  have  a  >  00(100?);  thus  u^  >  . 

From  (ii),  1  >  10  <*  •  By  (iii)  and  (iv),  o010a  >o01>  1  .  So 

u^  >  Ug  >  1  . 


Now  assume  u^_1  >  u^  >  1  .  By  (ii),  1  >  10u^  >  •  By 

(iii)  and  (iv),  a  0  (l  /  u^)  >  or  0  (l  /  >  or©  1  >  1  .  So, 

>  1  •  By  induction,  the  sequence  [u^]  is  bounded  and  monotone. 
Therefore,  since  there  are  a  finite  number  of  floating-point  representations 
between  a  and  1  ,  the  sequence  converges  to  a  limit  u  >  1  .  In  the 
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limit,  we  have 

u  =  Of©  (I0u)  . 

Following  the  techniques  o  f  Wilkinson  (19&3),  we  have 
u  =  (a  -  u"1  (i+c))(i+Tl) 
for  some  values  of  s  and  T|  satisfying 

|  c|  <  Pi_t  and  |H  |  <  a1_t  . 


u  =  or  -  u-1  +  6  , 


where  6  =  QtT\  -  u-i(c  +  T|  +  *T|)  • 
Therefore, 


u  =£[(<*  +  6)  + y(a  +6)  -  4  ] . 


From  Theorem  1  we  see  that 


u  -  u  =  0(6)  =  0(8  ) 


Since  u  >  1  ,  Theorem  4  provides  a  bound  on  the  relative  error  in  u  . 

We  would  like  to  remark  that  the  algorithm  (2)  is  nothing  more  than 
Gaussian  elimination  which  is  known  to  be  very  stable  for  positive  definite 
systems.  The  condition  number  of  the  matrix  B  is  easily  calculated  to  be 


•ond(B)  = 


I  a  I  +2  cos 


lo  I  -  2  cos 


a\  +  2 

lor  I  -  2 


Using  the  error  bound  given  in  Forsythe  and  Moler  (1967):  If 
By  =  d  and  (B+E)z  =  d  ,  then 


lly  -  51!  PI! 

-  <  cond(B)  -  , 

ll£ll  l|B|! 
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where  ||*||  denotes  the  spectral  norm.  If  E  is  due  to  roundoff  error 
in  representing  a  ,  then  ||E||  <  e  =  |ajp1-t  ,  and 


lly  -  ?ll 
Il5ll 


c 

| a  |  -  2  cos 


n 

n+1 


h .  Generalizations 

An  important  extension  of  Theorem  1  is  that  the  LU  decomposition 
will  converge  even  if  some  of  the  upper  left  elements  of  the  matrix  are 
changed.  If  a  tri-diagonal  matrix  contains  a  Toeplitz  sub-matrix,  then 
that  portion  of  the  LU  decomposition  converges.  Problems  of  this  sort 
occur,  for  example,  with  cubic  spline  interpolation  with  prescribed  deriv¬ 
atives  at  the  ends.  This  is  a  result  of  the  following. 

Theorem  5 :  If  oi>2  and  u1  =  Y  where  Y  has  any  value  except  0  , 
l/a  ,  or  u  ,  then  the  sequence  ui  =  a  -  l/u^  ^  ,  i=2, . . . ,  converges 
to  u+  where 

Of  +~V  a  -  4 

u  =  Jr.  i. 

+  2  » 

and  _____ 

»  -Vo2  -  4 


(A  similar  result  holds  for  0/  <  -2  .) 

Proof;  The  nonlinear  difference  equation,  u^^ 


explicitly  by  using  the  substitution  u.  =  - 

wi-l 

second-order  difference  equation.  For  a  >  0 


=  Q 1 - ,  can  be  solved 

ui-l 

■  to  produce  a  linear 
and  u^  =  Y  t  the  solution 


is; 
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ui  =  u+ 


1  +  ? 


'  v  i+1 

(t) 


where  ?  = 


V 


_  Vq,2  -  4  -  V  +  u. 


V  -  u 


.  Since  a  >  2  ,  the  positive  quantity 


(u  /u+)  is  less  than  unity.  Convergence  follows  immediately. 


The  results  we  have  given  for  scalars  can  also  be  generalized  to 
matrices . 

Theorem  6;  If  a  matrix  can  be  partitioned  as 

A  B 
BAB 

BAB 


a  = 


where  both  A  and  B  are  symmetric  and  positive  definite,  and  if  the 
eigenvalues  of  B_1A  are  greater  than  2  in  modulus ,  then  the  block 
Gaussian  elimination  of  d  converges . 

Proof :  Block  elimination  is  equivalent  to  constructing  the  sequence  of 

matrices  U  =  A  ,  Ui+1  =  A  -  BU^B  ,  i=l,2,...  .  But  A  =  PAPT  and 

T  -1 

B  =  PP  where  A  is  the  diagonal  matrix  of  eigenvalues  of  B  A  .  Define 


A  =  A  and  A^+1 


=  A  -  A’1.  Then  Ux  =  PA-,?1  and  if  =  PAiPT  then 


Ui+1  =  PApT  "  PPT(PAiPT)“1PPT 

-IT  T 

=  P[ A  -  A.  ]PA  =  PAi+1P 


The 


convergence  of  A^  (as  well  as  the  rate  of  convergence)  under  the 
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conditions  stated  follows  from  the  results  for  scalars  given  in  Theorems  1-3 
An  example  of  a  matrix  that  satisfies  the  required  conditions  for 
convergence  is  the  matrix  that  arises  from  the  five-point  finite  difference 
approximation  to  Laplace's  operator'  in  a  rectangle: 

n  -i  1 


-i 

4  -l 

•  • 

•  •  • 

•  •  • 

-1 

-1  4 

However,  this  method  does  not  appear  to  be  competitive  with  existing 
methods  for  this  particular  matrix. 

5 .  Conclusions 

Many  of  the  observations  which  lead  to  the  simplification  in  com¬ 
puting  the  LU  decomposition  for  tri -diagonal  Toepli tz  matrices  generalize 
to  Toeplitz  band  matrices.  Bauer  (1955)  states  that  the  Cholesky  decomp¬ 
osition  of  band  symmetric  matrices  converges  in  the  sense  that  each 
diagonal  of  the  triangular  matrix  converges.  We  know  of  no  rate-of- 
convergence  results  for  the  band  case. 

An  alternate  proof  of  Theorem  1  can  be  easily  constructed  by  con¬ 
sidering  the  analytical  solution  to  the  difference  equation  (l).  Bounds 


where 


A  = 


15 


on  k  similar  to  those  given  in  Theorems  2  and  3,  but  not  quite  as  close 
can  be  obtained  similarly. 


14 


I 

! 


Acknowledgement 


We  would  like  to  thank  Professors  Gene  H.  Golub  and  Gerald  Taylor 
for  criticizing  the  manuscript  and  for  several  stimulating  discussions. 
In  addition  we  thank  Professor  Golub  for  bringing  the  work  of  Parter  and 
Bauer  to  our  attention  and  for  suggesting  the  proof  technique  used  in 
Theorem  4 . 


15 


REFERENCES 


Bauer,  Friedrich  L.  (1955),  "Ein  direktes  Iterationsverfahren  zur 
Hurwitz-Zerlegung  eines  Polynoms,"  Archiv  der  Elektrischen  Uebertragung, 
Wiesbaden,  Vol.  9,  285-290.  Translated  from  the  German  by  M.  Morf. 

Forsythe,  G.  E.  and  Moler,  C.  B.  (1967),  Computer  Solution  of  Linear 
Algebraic  Equations,  Englewood  Cliffs ,  N.J.:  Prentice-Hall,  Inc. 

Parter,  Seymour  V.  (1962),  "An  Observation  on  the  Numerical  Solution  of 
Difference  Equations  and  a  Theorem  of  Szego,"  Num.  Math.  4,  295-295. 

Wilkinson,  J.  H.  (19^5 ),  Rounding  Errors  in  Algebraic  Processes. 
Englewood  Cliffs,  N.J.:  Prentice-Hall,  Inc. 


16 


